Jupyter Notebook for Julia¶
A. Install Julia in your Google Colab.¶
If Julia is already isntalled in your Google Colab, go to B. If Julia and Jupyter Notebook are installed in your computer, go to D.
"Run" means
- Put the cursor on the cel and shift+enter on Jupyter Notebook
- Tap the circle near the upper left corner of the cell on Google Colab
Run the following cell.¶
!curl -fsSL https://install.julialang.org | sh -s -- --yes
!julia -v
using Pkg
Pkg.add("LaTeXStrings")
Pkg.add("Images")
Pkg.add("Plots")
Pkg.add("Colors")
Pkg.add("Downloads")
Pkg.add("FileIO")
Pkg.add("LowRankApprox")
Pkg.add("Wavelets")
Pkg.add("DSP")
Pkg.add("Dates")
Pkg.add("CSV")
Pkg.add("DataFrames")
Pkg.add("Statistics")
using Random
using LinearAlgebra
using LaTeXStrings
using Images
using Plots
using Colors
using Downloads
using FileIO
using LowRankApprox
using Wavelets
using DSP
using Dates
using CSV
using DataFrames
using Statistics
# space
Topics¶
- The Collatz conjecture
- RGB images
- RGB movies
- Singular value decomposition and low rank approximation
- Wavelet decomposition
- Public data and visualization
- Central limit theorem for the binomial distribution $B(n,p)$
- Taylor series of $\sin{x}$
- Fourier series of a step function and a sawtooth function
- Weyl's equidistribution theorem
1. The Collatz conjecture¶
Generate a sequence $\{x[n]\}_{n=1}^\infty$ by the following process:
- Given an arbitrary positive integer $m$, set $x[1]:=m$.
- For $n=2,3,4,\dotsc$, if $x[n-1]:=1$, stop the process, otherwise
\begin{aligned} x[n] & :=\dfrac{x[n-1]}{2}\quad (x[n-1]\ \text{is even}), \\ x[n] & := 3x[n-1]+1\quad (x[n-1]\ \text{is odd}). \end{aligned}
In 1937 Lothar Collatz introduce the following conjecture:
For any $m=1,2,3,\dotsc$, there exists $n_m=1,2,3,\dotsc$ such that $x[n_m]=1$.¶
See this paper for instance. Here are some examples:
$1, \quad 2\rightarrow1, \quad 3\rightarrow10\rightarrow5\rightarrow16\rightarrow8\rightarrow4\rightarrow2\rightarrow1, \quad 4\rightarrow2\rightarrow1, \quad\dotsc$
Barina verified the Collatz conjecture is true for all $m\leqq2^{68}$ using workstations. See this paper.
We examine this conjecture for $m=1,2,3,\dotsc,500$.
Run the following cell to create Collatz sequences for $x[1]=m=1,\dotsc,500$, and observe the behavior of each sequence.¶
# upper bound of m
M=500;
x=ones(Int64,150);
animcollatz = @animate for m=1:M
x[1]=m
for n=2:150
if x[n-1]==1
x[n]=x[n]
else
if mod(x[n-1],2)==0
x[n]=Int64.(x[n-1]/2)
else x[n]=3*x[n-1]+1
end
end
end
scatter(x,
ylim=(1,2000),
title="Collatz sequence x[n] with x[1]=$m",
xlabel=L"n",
ylabel=L"x[n]",
xticks=([1,20,40,60,80,100,120,140]),
legend=false,
clear=true)
end
gif(animcollatz, "collatz_animation.gif", fps = 5)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/collatz_animation.gif
Run the following cell to plot the pairs of $(m,n_m)$ for $m=1,\dotsc,500$ in a scatter graph.¶
function collatz(n)
n % 2 == 0 ? n ÷ 2 : 3n + 1
end
function stopping_time(m)
x, steps = m, 0
while x > 1
x = collatz(x)
steps += 1
end
return steps
end
N = 1:500
steps = [stopping_time(m) for m in N]
scatter(N, steps,
legend=false,
xlabel=L"\mathrm{Initial}\ \mathrm{value}\ \ m",
ylabel=L"\mathrm{Steps}\ \ n_m\ \ \mathrm{to}\ \ 1",
title="Collatz stopping steps", markersize=2)
2. RGB images¶
The data of an RGB image is a triplet of matrices of the same size. This is essentially a triplet of grayscale images of the same size. These three grayscale images, colored red, green, and blue, are then overlaid to form an RGB image.
Run the following cell to create three $m{\times}n=10\times15$ matrices $R,G,B$ by random numbers.¶
Random.seed!(1)
m, n = 10, 15 # size of matrices
Z = zeros(m, n); # zero matrix
R = rand(m, n);
G = rand(m, n);
B = rand(m, n) # See below
10×15 Matrix{Float64}:
0.670414 0.643352 0.376138 0.983634 … 0.899505 0.380349 0.221796
0.691898 0.100056 0.218766 0.40892 0.424926 0.958143 0.950771
0.192718 0.164962 0.60248 0.0424498 0.237027 0.610986 0.839523
0.516169 0.959669 0.753166 0.705895 0.189413 0.619391 0.736151
0.0475605 0.72755 0.852577 0.0647826 0.049791 0.849654 0.777125
0.245685 0.436295 0.88632 0.792918 … 0.587347 0.515305 0.0444534
0.979139 0.396001 0.0642346 0.080433 0.171102 0.0882561 0.0745833
0.0294498 0.419402 0.890999 0.998768 0.265295 0.811548 0.626937
0.224011 0.525583 0.988467 0.0620525 0.480608 0.0926662 0.601864
0.80577 0.21062 0.292272 0.199147 0.582327 0.494502 0.440245
Run the following cell to colorize $R,G,B$, and show red, green, blue and RGB images.¶
imgR = colorview(RGB, R, Z, Z)
imgG = colorview(RGB, Z, G, Z)
imgB = colorview(RGB, Z, Z, B)
imgRGB = colorview(RGB, R, G, B)
p1 = plot(imgR, grid=false, axis=false, title="Red")
p2 = plot(imgG, grid=false, axis=false, title="Green")
p3 = plot(imgB, grid=false, axis=false, title="Blue")
p4 = plot(imgRGB, grid=false, axis=false, title="RGB")
plot(p1, p2, p3, p4, layout=(2,2), size=(600,450))
Random.seed!(1)
m, n = 10, 15 # size of matrices
K = 100 # steps, length of the sequence
Z = zeros(m, n);
RR = rand(m, n, K);
GG = rand(m, n, K);
BB = rand(m, n, K);
Run the following cell to see the movie.¶
animrgb = @animate for k=1:K
movieR = colorview(RGB, RR[:,:,k], Z, Z)
movieG = colorview(RGB, Z, GG[:,:,k], Z)
movieB = colorview(RGB, Z, Z, BB[:,:,k])
movieRGB = colorview(RGB, RR[:,:,k], GG[:,:,k], BB[:,:,k])
p1 = plot(movieR, grid=false, axis=false, title="Red")
p2 = plot(movieG, grid=false, axis=false, title="Green")
p3 = plot(movieB, grid=false, axis=false, title="Blue")
p4 = plot(movieRGB, grid=false, axis=false, title="RGB")
plot(p1, p2, p3, p4, layout=(2,2), size=(600,450))
end
gif(animrgb, "rgb_animation.gif", fps = 3)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/rgb_animation.gif
4. Singular value decomposition and low rank approximation¶
In general an $m{\times}n$ matrix $A=\begin{bmatrix}\vec{a}_1 & \dotsb & \vec{a}_n \end{bmatrix}$ has a specific nonnegative integer which is called as the rank of $A$ defined by
$$r := \operatorname{dim}\bigl( \{ A\vec{x} = x_1\vec{a}_1+\dotsb+x_n\vec{a}_n \ \vert \ \vec{x} \in \mathbb{R}^n \} \bigr).$$
It follows that $1 \leqq r \leqq \min\{m,n\}$ unless all the entries of $A$ are zero. It is known that there exist
$$\sigma_1 \geqq \dotsb \geqq \sigma_r>0, \quad \vec{u}_1,\dotsc,\vec{u}_r \in \mathbb{R}^m, \quad \vec{v}_1,\dotsc,\vec{v}_r \in \mathbb{R}^n$$
such that we have the singular value decomposition:
$$A=\sum_{j=1}^r\sigma_j\vec{u}_j\vec{v}_j^T.$$
Set
$$A_k:=\sum_{j=1}^k\sigma_j\vec{u}_j\vec{v}_j^T, \quad k=1,\dotsc,r.$$
It is known that the rank of $A_k$ is $k$ and closest to $A$ among all the $m{\times}n$ matrices of rank $k$, that is, $A_k$ is the best approximation of $A$ among all the $m{\times}n$ matrices of rank $k$.
We compute a $432\times768\times3$ RGB image of hibiscus. All the R, G, B matrices are of rank $432$.
Run the following cell to compute the SVD of the RGB-elements of the picture.¶
url = "https://raw.githubusercontent.com/fiomfd/hands-on/refs/heads/main/data/hibiscus.jpg"
Downloads.download(url, "hibiscus.jpg")
I=load("./hibiscus.jpg");
X=imresize(I, ratio=1/5);
(p,q)=size(X);
AAA=channelview(X);
RRR=Array{Float64}(AAA[1,:,:]);
GGG=Array{Float64}(AAA[2,:,:]);
BBB=Array{Float64}(AAA[3,:,:]);
RU, RS, RV=psvd(RRR);
GU, GS, GV=psvd(GGG);
BU, BS, BV=psvd(BBB);
rank=60;
DR=zeros(p,q,rank);
DG=zeros(p,q,rank);
DB=zeros(p,q,rank);
for r=1:rank
DR[:,:,r]=sum(RS[n]*RU[1:p,n]*(RV[1:q,n])' for n=1:r);
DG[:,:,r]=sum(GS[n]*GU[1:p,n]*(GV[1:q,n])' for n=1:r);
DB[:,:,r]=sum(BS[n]*BU[1:p,n]*(BV[1:q,n])' for n=1:r);
end
Run the following cell to see the movie of the best low rank approximation for $r=1,\dotsc,60$.¶
animsvd = @animate for k=1:rank
W=zeros(3,p,q,rank);
for r=1:rank
W[1,:,:,r]=DR[:,:,r];
W[2,:,:,r]=DG[:,:,r];
W[3,:,:,r]=DB[:,:,r];
end
SVD1=plot(colorview(RGB,W[:,:,:,k]),
title="Rank-$k approximation",
titlefontsize=12,
xaxis=false,
xticks=false,
yaxis=false,
yticks=false,
grid=false);
SVD2=plot(colorview(RGB,X),
title="Original RGB image (rank 432)",
titlefontsize=12,
xaxis=false,
xticks=false,
yaxis=false,
yticks=false,
grid=false);
plot(SVD1,SVD2, layout=(1,2), size=(650,220))
end
gif(animsvd, "svd_animation.gif", fps = 3)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/svd_animation.gif
5. Wavelet decomposition¶
The discrete wavelet is a pair of orthonomal two vectors $\vec{u}$ and $\vec{v}$ satisfying some condition, and the most typical example of the discrete wavelets is the Haar wavelet. They are very simple vectors:
$$\vec{u} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ 1 \\ 0 \\ \vdots \\ 0 \end{bmatrix}, \quad \vec{v} = \frac{1}{\sqrt{2}} \begin{bmatrix} 1 \\ -1 \\ 0 \\ \vdots \\ 0 \end{bmatrix} \in\mathbb{R}^{N},$$
Roughly speaking the discrete wavelet transform by the Haar wavelet is to take averages of neighboring elements. The elements of a vector are transformed to the average of neighboring two elements, and the elements of a matrix are converted to the average of neighboring 2 by 2 elements. The following examples show the outcomes of repeating the discere wavelet tranform twice:
$$\vec{a}_0 = \begin{bmatrix} 1 \\ 3 \\ 5 \\ 7 \end{bmatrix} \mapsto \vec{a}_1 = \begin{bmatrix} 2 \\ 2 \\ 6 \\ 6 \end{bmatrix} \mapsto \vec{a}_2 = \begin{bmatrix} 4 \\ 4 \\ 4 \\ 4 \end{bmatrix}$$
$$A_0 = \begin{bmatrix} 0 & 2 & 4 & 6 \\ 8 & 10 & 12 & 14 \\ 16 & 18 & 20 & 22 \\ 24 & 26 & 28 & 30 \end{bmatrix} \mapsto A_1 = \begin{bmatrix} 5 & 5 & 9 & 9 \\ 5 & 5 & 9 & 9 \\ 21 & 21 & 25 & 25 \\ 21 & 21 & 25 & 25 \end{bmatrix} \mapsto A_2 = \begin{bmatrix} 15 & 15 & 15 & 15 \\ 15 & 15 & 15 & 15 \\ 15 & 15 & 15 & 15 \\ 15 & 15 & 15 & 15 \end{bmatrix}.$$
The output such as $\vec{a}_\ell$ and $A_\ell$ of $\ell$ times operations is said to be the approximation part of level $\ell$, and the remainder term such as $\vec{a}_0-\vec{a}_\ell$ and $A_0-A_\ell$ is called the detail part of level $\ell$. In terms of Fourier analysis, the approximation part is the low frequency part, and the detail part is the high frequency part and is similar to the edge detection.
We observe the wavelet decomposition of the $384\times512$ RGB image of fried noodle using the Haar wavelet. We can get the wavelet decomposition up to level $7$ since $384=2^7\times3$ and $512=2^9$.
Run the following cell to implement the wavelet decomposition.¶
url = "https://raw.githubusercontent.com/fiomfd/hands-on/refs/heads/main/data/CityU.jpg"
Downloads.download(url, "CityU.jpg")
I8=load("./CityU.jpg");
X8=imresize(I8, ratio=1/8);
(p8,q8)=size(X8);
A8=channelview(X8);
R8=Array{Float64}(A8[1,:,:]);
G8=Array{Float64}(A8[2,:,:]);
B8=Array{Float64}(A8[3,:,:]);
# Maximum level
L=7;
# Decomposition Filter
XR=zeros(p8,q8,L);
XG=zeros(p8,q8,L);
XB=zeros(p8,q8,L);
for l=1:L
XR[:,:,l]=dwt(R8, wavelet(WT.haar), l);
XG[:,:,l]=dwt(G8, wavelet(WT.haar), l);
XB[:,:,l]=dwt(B8, wavelet(WT.haar), l);
end
# Splitting
XRapprox=zeros(p8,q8,L);
XGapprox=zeros(p8,q8,L);
XBapprox=zeros(p8,q8,L);
for l=1:L
XRapprox[1:Int(p8/2^l),1:Int(q8/2^l),l]=XR[1:Int(p8/2^l),1:Int(q8/2^l),l];
XGapprox[1:Int(p8/2^l),1:Int(q8/2^l),l]=XG[1:Int(p8/2^l),1:Int(q8/2^l),l];
XBapprox[1:Int(p8/2^l),1:Int(q8/2^l),l]=XB[1:Int(p8/2^l),1:Int(q8/2^l),l];
end
XRdetail=XR-XRapprox;
XGdetail=XG-XGapprox;
XBdetail=XB-XBapprox;
# Composition Filter
YRapprox=zeros(p8,q8,L);
YGapprox=zeros(p8,q8,L);
YBapprox=zeros(p8,q8,L);
YRdetail=zeros(p8,q8,L);
YGdetail=zeros(p8,q8,L);
YBdetail=zeros(p8,q8,L);
for l=1:L
YRapprox[:,:,l]=idwt(XRapprox[:,:,l], wavelet(WT.haar), l);
YGapprox[:,:,l]=idwt(XGapprox[:,:,l], wavelet(WT.haar), l);
YBapprox[:,:,l]=idwt(XBapprox[:,:,l], wavelet(WT.haar), l);
YRdetail[:,:,l]=idwt(XRdetail[:,:,l], wavelet(WT.haar), l);
YGdetail[:,:,l]=idwt(XGdetail[:,:,l], wavelet(WT.haar), l);
YBdetail[:,:,l]=idwt(XBdetail[:,:,l], wavelet(WT.haar), l);
end
# RGB
Wapprox=zeros(3,p8,q8,L+1);
Wdetail=zeros(3,p8,q8,L+1);
Wapprox[:,:,:,1]=A8;
for l=2:L+1
Wapprox[1,:,:,l]=YRapprox[:,:,l-1];
Wapprox[2,:,:,l]=YGapprox[:,:,l-1];
Wapprox[3,:,:,l]=YBapprox[:,:,l-1];
Wdetail[1,:,:,l]=YRdetail[:,:,l-1];
Wdetail[2,:,:,l]=YGdetail[:,:,l-1];
Wdetail[3,:,:,l]=YBdetail[:,:,l-1];
end
Run the previous cell to see a movie of the wavelet decomposition.¶
animhaar = @animate for l=0:L
P1=plot(colorview(RGB,A8[:,:,:]),
title="Original RGB image",
titlefontsize=12,
xaxis=false,
xticks=false,
yaxis=false,
yticks=false,
grid=false);
P2=plot(colorview(RGB,Wapprox[:,:,:,l+1]),
title="Level $l Approximation",
titlefontsize=12,
xaxis=false,
xticks=false,
yaxis=false,
yticks=false,
grid=false);
P3=plot(colorview(RGB,Wdetail[:,:,:,l+1]),
title="Level $l Detail",
titlefontsize=12,
xaxis=false,
xticks=false,
yaxis=false,
yticks=false,
grid=false);
plot(P1,P2,P3,size=(900,240),layout=(1,3))
end
gif(animhaar, "svd_animation.gif", fps = 0.2)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/svd_animation.gif
6. Public data and visualization¶
Let $\vec{x}=\bigl[x[1],\dotsc,x[N]\bigl]^T\in\mathbb{R}^N$ be an $N$-dimensional real column vector. This can describes data such as $x[i]$ means the highest temperature of the $i$-th day for $i=1,\dotsc,N$. One can exploit some information such as the average $m$ and the variance $\sigma^2$:
$$m:=\frac{1}{N}\sum_{i=1}^Nx[i], \quad \sigma^2:=\frac{1}{N}\sum_{i=1}^N(x[i]-m)^2.$$
In what follows Julia downloads open data, and visualizes information taking averages, maximums and minimums. More concretely Julia downloads
- the daily maximum tenperature
- the daily minimum tenpereture
- the daily mean tenperature of 1884-present at the Hong Kong Observatory. The starting dates of the theree data are not same, and so we make use of data of 140 years of 1885-2024. The Hong Kong Observatory dataset has no daily temperature records from 1 Jan 1940 to 31 Dec 1946, due to the disruption of WWII and the Japanese occupation. This 7-year gap means that analyses requiring continuous daily data over multiple years will be affected if these years are included. To ensure consistency, we restrict the analysis period to 1 Jan 1947 – 31 Dec 2024.
Run the following cell to donwload the public data and show the table.¶
Downloads.download("https://data.weather.gov.hk/weatherAPI/opendata/opendata.php?dataType=CLMMAXT&rformat=csv&station=HKO","HK_Maximun_Temperature.csv");
Downloads.download("https://data.weather.gov.hk/weatherAPI/opendata/opendata.php?dataType=CLMMINT&rformat=csv&station=HKO","HK_Minimum_Temperature.csv");
Downloads.download("https://data.weather.gov.hk/weatherAPI/opendata/opendata.php?dataType=CLMTEMP&rformat=csv&station=HKO","HK_Mean_Temperature.csv");
Amax = CSV.read("HK_Maximun_Temperature.csv", DataFrame,
header=["Year", "Month", "Day", "Maximum", "C"],
skipto=370,
footerskip=3,
types=Dict("Maximum" => Float64),
drop=["C"],
missingstring="***",
normalizenames=true);
delete!(Amax, 48579:nrow(Amax));
Amean = CSV.read("HK_Mean_Temperature.csv", DataFrame,
header=["Year", "Month", "Day", "Average", "C"],
skipto=310,
footerskip=3,
types=Dict("Average" => Float64),
drop=["C"],
missingstring="***",
normalizenames=true);
delete!(Amean, 48579:nrow(Amean));
Amin = CSV.read("HK_Minimum_Temperature.csv", DataFrame,
header=["Year", "Month", "Day", "Minimum", "C"],
skipto=370,
footerskip=3,
types=Dict("Minimum" => Float64),
drop=["C"],
missingstring="***",
normalizenames=true);
delete!(Amin, 48579:nrow(Amin));
A = outerjoin(Amax, Amean, Amin,on = [:Year, :Month, :Day])
| Row | Year | Month | Day | Maximum | Average | Minimum |
|---|---|---|---|---|---|---|
| Int64 | Int64 | Int64 | Float64? | Float64? | Float64? | |
| 1 | 1885 | 1 | 1 | 17.8 | 15.3 | 13.9 |
| 2 | 1885 | 1 | 2 | 20.5 | 16.7 | 14.2 |
| 3 | 1885 | 1 | 3 | 16.6 | 15.4 | 14.8 |
| 4 | 1885 | 1 | 4 | 18.6 | 16.1 | 14.6 |
| 5 | 1885 | 1 | 5 | 20.8 | 17.5 | 15.1 |
| 6 | 1885 | 1 | 6 | 21.1 | 18.3 | 16.8 |
| 7 | 1885 | 1 | 7 | 19.9 | 17.9 | 15.4 |
| 8 | 1885 | 1 | 8 | 17.3 | 15.6 | 14.2 |
| 9 | 1885 | 1 | 9 | 16.9 | 15.7 | 14.7 |
| 10 | 1885 | 1 | 10 | 18.1 | 16.1 | 15.1 |
| 11 | 1885 | 1 | 11 | 19.9 | 17.2 | 15.3 |
| 12 | 1885 | 1 | 12 | 17.9 | 16.0 | 14.9 |
| 13 | 1885 | 1 | 13 | 17.1 | 14.4 | 11.9 |
| ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ | ⋮ |
| 48567 | 2024 | 12 | 20 | 17.7 | 14.9 | 11.9 |
| 48568 | 2024 | 12 | 21 | 20.2 | 16.9 | 13.9 |
| 48569 | 2024 | 12 | 22 | 18.0 | 15.8 | 13.5 |
| 48570 | 2024 | 12 | 23 | 17.5 | 16.5 | 15.1 |
| 48571 | 2024 | 12 | 24 | 19.1 | 17.4 | 15.6 |
| 48572 | 2024 | 12 | 25 | 20.6 | 18.5 | 16.6 |
| 48573 | 2024 | 12 | 26 | 22.9 | 20.1 | 18.0 |
| 48574 | 2024 | 12 | 27 | 20.9 | 19.2 | 18.1 |
| 48575 | 2024 | 12 | 28 | 18.8 | 16.9 | 15.1 |
| 48576 | 2024 | 12 | 29 | 17.4 | 15.4 | 13.3 |
| 48577 | 2024 | 12 | 30 | 20.4 | 17.7 | 14.3 |
| 48578 | 2024 | 12 | 31 | 22.6 | 19.8 | 17.6 |
Run the following cell to draw the line graphs of Annual maximum, average, and minimum temperatures.¶
Over the long term, temperatures are rising slowly.
A_1947_2024 = filter(:Year => y -> 1947 <= y <= 2024, A);
Y = combine(groupby(A_1947_2024, :Year),
:Maximum => maximum => :Ymax,
:Average => mean => :Yavg,
:Minimum => minimum => :Ymin
)
plot(Matrix(Y)[:,2:4],
size=(500,350),
title="Annual max, ave, and min temperatures at HK Observatory",
titlefontsize=9,
titlefontcolor="blue",
xlim=(-6,84),
xticks=([4,14,24,34,44,54,64,74],
[1950 1960 1970 1980 1990 2000 2010 2020]),
ylim=(0,40),
grid=false,
xlabel="year",
ylabel="C",
legend=false)
Run the following cell to create a movie showing the monthy maximum, average, and minmum temperatures in each year between 1947 and 2024.¶
animhko = @animate for k=1947:2024
M = combine(groupby(A_1947_2024, [:Year, :Month]),
:Maximum => maximum => :Mmax,
:Average => mean => :Mavg,
:Minimum => minimum => :Mmin)
Myy = sort(M[M.Year .== k, :], :Month)
plot(Matrix(Myy[:,3:5]),
size=(500,350),
title="Monthly max, ave, and min temperatures at HK Observatory in $k",
titlefontsize=8,
titlefontcolor="blue",
xlim=(0,13),
xticks=([1,2,3,4,5,6,7,8,9,10,11,12]),
ylim=(0,40),
grid=false,
xlabel="month",
ylabel="C",
legend=false)
end
gif(animhko, "HKO_animation.gif", fps = 2)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/HKO_animation.gif
7. Central limit theorem for the binomial distribution $B(n,p)$¶
Let $n$ be a positive integer, and let $p$ be a constant satisfying $0 < p < 1$. Denote by $B(n,p)$ the the binomial distribution of $n$ iid experiments $X_1,\dotsc,X_n$ with probability $p$, that is, a discrete probability distribution on $\{0,1,\dotsc,n\}$ such that
$$\operatorname{P}(S_n=k)=\frac{n!}{k!(n-k)!}p^k(1-p)^{n-k}, \quad k=0,1,\dotsc,n,$$
where $S_n:=X_1+\dotsb+X_n$.
It is well-known that for any $a,b\in\mathbb{R}$ with $a<b$,
$$\operatorname{P} \left( a \leqq \frac{S_n-pn}{\sqrt{np(1-p)}} \leqq b \right) \rightarrow \frac{1}{\sqrt{2\pi}} \int_a^be^{-x^2/2}dx \quad (n \rightarrow \infty),$$
which is a typical example of the central limit theorem.
Choose $p\in(0,1)$ and run the following cell to create a movie of the central limit theorem¶
# Choose 0<p<1
p=0.3;
x1 = range(-3, 3, length = 1201);
f1 = exp.(-x1.^2/2)/sqrt(2*pi);
animclt = @animate for m=1:150
n1=10*m
z2=zeros(n1+2*max(100,n1));
for l=0:n1
z2[l+max(100,n1)+1]=Float64(sqrt(n1*p*(1-p))*factorial(BigInt(n1))*p^l*(1-p)^(n1-l)/factorial(BigInt(l))/factorial(BigInt(n1-l)));
end
z1 = zeros(length(x1));
c1 = Int64((length(x1)-1)/2);
for k=1:c1
z1[c1+k]=z2[max(100,n1)+2+Int64(floor(n1*p+(k-1)/200*sqrt(n1*p*(1-p))))];
z1[c1+1-k]=z2[max(100,n1)+1+Int64(floor(n1*p-(k-1)/200*sqrt(n1*p*(1-p))))];
end
plot(f1,
title=(L"B(n,%$p) \rightarrow N(0,1), n=%$n1"),
size=(500,350),
grid=false,
ylim=(0,0.45),
xticks=([201 401 601 801 1001;],[-2,-1,0,1,2]),
yticks=([0 1/sqrt(2*pi);],[0,"1/√2π"]),
xlabel=(L"x"),
ylabel=(L"y"),
edgecolor=:false,
legend=:false,
lw=2,
color=:blue)
plot!(z1, seriestype=:shape, color=:yellow, alpha=0.6)
end
gif(animclt, "CLT_animation.gif", fps = 10)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/CLT_animation.gif
8. Taylor series of $\sin{x}$¶
It is well-known that $\sin{x}$ ($x\in\mathbb{R}$) equals a power series of $s$:
$$\sin{x}=\sum_{k=0}^\infty\frac{(-1)^kx^{2k+1}}{(2k+1)!}, \quad x\in\mathbb{R}.$$
The right hand side is said to be the Taylor series of $\sin{x}$ or the Taylor expansion of $\sin{x}$. More precisely, for any $R>0$,
$$\max_{x\in[-R,R]} \left\lvert \sin{x} - \sum_{k=0}^K \frac{(-1)^kx^{2k+1}}{(2k+1)!} \right\rvert \rightarrow 0 \quad (K\rightarrow\infty).$$
We now observe this convergence. Set
$$S_K(x) := \sum_{k=0}^K \frac{(-1)^kx^{2k+1}}{(2k+1)!}.$$
Run the following cell to create a movie.¶
x = range(-3*pi, 3*pi, length = 301);
f3=sin.(x);
KK=14;
s3=zeros(KK+1,301);
for k=0:KK
for l=1:301
s3[k+1,l]=(-1)^k*x[l]^(2*k+1)/factorial(big(2*k+1));
end
end
S3=zeros(KK+1,301);
S3[1,:]=s3[1,:];
for l=1:301
for k=2:KK+1
S3[k,l]=S3[k-1,l]+s3[k,l];
end
end
animtaylor = @animate for K=0:KK
h3=S3[K+1,:];
plot([f3,h3],grid=false,linewidth=2,ylim=(-1.2,1.2),
size=(500,350),
title=(L"\sin(x), S_K(x), \quad K=%$K"),
xticks = ([0 50 100 151 201 251 301;], ["-3π","-2π","-π","0","π","2π","3π"]),
xlabel=L"x",
label=[L"\sin(x)" L"S_K(x)"],
legend=:topleft,legendfont=font(10))
end
gif(animtaylor, "Taylor_animation.gif", fps = 1)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/Taylor_animation.gif
9. Fourier series of a step function and a sawtooth function¶
Let $g(x)$ and $h(x)$ be a $1$-periodic step function and a $1$ periodic sawtooth function, that is,
$$\begin{aligned} g(x) & := \begin{cases} 1, &\ x\in[0,1/2], \\ 0, &\ x\in(1/2,1), \end{cases} \qquad g(x):=g(x-n),\ x\in[n,n+1),\ n\in\mathbb{Z}, \\ h(x) & :=x, \qquad h(x):=h(x-n),\ x\in[n,n+1),\ n\in\mathbb{Z}. \end{aligned}$$
Their Fourier series are
$$\begin{aligned} g(x) & \sim \frac{1}{2} + \sum_{k=1}^\infty \frac{2}{(2k-1)\pi} \sin\bigl(2(2k-1)\pi x\bigr), \\ h(x) & \sim \frac{1}{2} - \sum_{n=1}^\infty \frac{1}{n\pi} \sin(2\pi nx), \end{aligned}$$
Set
$$\begin{aligned} S_{2K-1}[g](x) & = \frac{1}{2} + \sum_{k=1}^K \frac{2}{(2k-1)\pi} \sin\bigl(2(2k-1)\pi x\bigr) \quad (K=1,2,3.\dotsc), \\ S_{-1}[g](x) & :=\frac{1}{2}, \\ S_N[h](x) & = \frac{1}{2} - \sum_{n=1}^N \frac{1}{n\pi} \sin(2\pi nx) \qquad (N=1,2,3,\dotsc), \\ S_0[h](x) & := \frac{1}{2}. \end{aligned}$$
On one hand, at the discontinuous points
\begin{aligned} S_{2K-1}[g](k) = \frac{1}{2} & \rightarrow \frac{1}{2} = \frac{g(k-0)+g(k+0)}{2} \quad (K\rightarrow\infty), \\ S_{2K-1}[g](k+1/2) = \frac{1}{2} & \rightarrow \frac{1}{2} = \frac{g(k+1/2-0)+g(k+1/2+0)}{2} \quad (K\rightarrow\infty), \\ S_N[h](k) = \frac{1}{2} & \rightarrow \frac{1}{2} = \frac{h(k-0)+h(k+0)}{2} \quad (N\rightarrow\infty) \end{aligned}
for all $k\in\mathbb{Z}$. On the other hand, we can prove that for any small $\delta>0$
\begin{aligned} \max_{x\in[\delta,1/2-\delta]\cup[1/2+\delta,1-\delta]} \lvert{S_{2K-1}[g](x)-g(x)}\rvert & \rightarrow 0 \quad (K\rightarrow\infty), \\ \max_{x\in[\delta,1-\delta]} \lvert{S_N[h](x)-h(x)}\rvert & \rightarrow 0 \quad (N\rightarrow\infty). \end{aligned}
The partial sums oscillate violately as $N$ increases near the discontinuous points. This is called the Gibbs phenomenon.
Run the following cell to create $g(x)$, $h(x)$, $S_{2K-1}[g](x)$ and $S_N[h](x)$.¶
# Set Maximun number of terms M2
M2=50;
# g
g = zeros(481)
for l=41:241
g[l] = 1;
end
for l=441:481
g[l]=1;
end
# S_{2K-1}[g]
x2 = range(-0.1, 1.1, length = 481);
S2=ones(M2+1,481)/2;
l=10
for k=2:M2+1
for l=1:481
S2[k,l]=S2[k-1,l]+2*sin(2*(2*k-3)*pi*x2[l])/(2*k-3)/pi;
end
end
#h
h = zeros(481);
for l=1:481
h[l]=x2[l]-Float64(floor(x2[l]));
end
# S_N[h]
S3=ones(M2+1,481)/2;
for k=2:M2+1
for l=1:481
S3[k,l]=S3[k-1,l]-sin(2*(k-1)*pi*x2[l])/(k-1)/pi;
end
end
Run the following cell to create a movie.¶
animgibbs = @animate for K2=1:M2+1
N2=K2-1
G1=plot(x2,S2[K2,:],
grid=false,
linewidth=2,
xlim=(-0.15,1.15),
ylim=(-1/8,9/8),
title="Step Function and its Fourier Series, K=$N2",
titlefontsize=9,
xticks = ([0 1/2 1;], [0,1/2,1]),
yticks = ([0 1/2 1;], [0,1/2,1]),
xlabel=L"x",
label=L"S_{2K-1}[g](x)",
legend=:topright,
legendfont=font(10))
plot!(x2,g,
linewidth=2,
label=L"g(x)")
G2=plot(x2,S3[K2,:],
grid=false,
linewidth=2,
xlim=(-0.15,1.15),
ylim=(-1/8,9/8),
title="Sawtooth Function and its Fourier Series, N=$N2",
titlefontsize=9,
xticks = ([0 1/2 1;], [0,1/2,1]),
yticks = ([0 1/2 1;], [0,1/2,1]),
xlabel=L"x",
label=L"S_N[h](x)",
legend=:top,
legendfont=font(10))
plot!(x2,h,
linewidth=2,
label=L"h(x)")
plot(G1,G2, layout=(1,2), size=(700,300))
end
gif(animgibbs, "Gibbs_animation.gif", fps = 5)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/Gibbs_animation.gif
10. Weyl's equidistribution theorem¶
Let $\lfloor{x}\rfloor$ be the floor function for $x\in\mathbb{R}$, that is, $$\lfloor{x}\rfloor:=\max\{m \in \mathbb{Z} : m \leqq x\},$$ and set
$$\langle{x}\rangle:=x-\lfloor{x}\rfloor.$$
Then $0\leqq\langle{x}\rangle<1$ and $x\equiv\langle{x}\rangle$ mod $1$.
We say that a sequence $\{a_n\}_{n=1}^\infty \subset [0,1)$ is equidistributed in $[0,1)$ if
$$\frac{\sharp\{n=1,\dotsc,N : a<a_n<b\}}{N} \rightarrow b-a \quad (N \rightarrow \infty)$$
Weyl's equidistribution theorem: For any $\gamma\in\mathbb{R}\setminus\mathbb{Q}$, $\{\langle{n\gamma}\rangle\}_{n=1}^\infty$ is equidistributed in $[0,1)$.
We observe Weyl's equidistribution theorem using a bijection
$$[0,1)\ni\langle{x}\rangle \mapsto e^{2\pi{i}\langle{x}\rangle} = e^{2\pi{i}x} = \cos(2\pi{x})+i\sin(2\pi{x}) \in \mathbb{S}^1:=\{z\in\mathbb{C} : \lvert{z}\rvert=1\}.$$
Choose $\gamma^2$ (positive integers, fractions and decimals), and run the following cell to create a movie.¶
# Choose g2=gamma^2
g2=10;
# Upper bound of N
Nmax=250
# Creating sequence
t10 = range(0, 1, length = 101);
c101=cos.(2*pi*t10);
c102=sin.(2*pi*t10);
γ=sqrt(g2);
x10=zeros(Nmax);
y10=zeros(Nmax);
for n=1:Nmax
x10[n]=cos(2*pi*γ*n);
y10[n]=sin(2*pi*γ*n);
end
animweyl = @animate for N10=1:Nmax
scatter(x10[1:N10],y10[1:N10],
title="Weyl's Equidistribution Theorem",
grid=false,
xlim=(-1.1,1.1),
ylim=(-1.1,1.1),
aspect_ratio=1.0,
xticks=false,
yticks=false,
label=false,
xaxis=false,
yaxis=false,
annotations=[(0,0.3,L"γ=\sqrt{%$g2}",20),(0,0,L"e^{2πiγn}, n=1,...,%$N10",20)],
annotationguide=:auto)
plot!(c101,c102,label=false,linecolor=:magenta)
end
gif(animweyl, "Weyl_animation.gif", fps = 3)
[ Info: Saved animation to /home/hc/Dropbox/ATCM/ATCM2026/Julia/Weyl_animation.gif